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The ground-state correlation functions of a one-dimensional homogeneous Bose system described 
by the Lieb-Liniger Hamiltonian are investigated by using exact quantum Monte Carlo techniques. 
This article is an extension of a previous study published in Phys. Rev. A 68, 031602 (2003). New 
results on the local three-body correlator as a function of the interaction strength are included and 
compared with the measured value from three-body loss experiments. We also carry out a thorough 
study of the short- and long-range behavior of the one-body density matrix. 
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I. INTRODUCTION 

Recent progress achieved in techniques of confining Bose condensates has lead to experimental realizations of 
quasi-onc- dimensional (ID) systcmsQ S S 0, S The quasi- ID regime is reached in highly anisotropic traps, 
where the axial motion of the atoms is weakly confined while the radial motion is frozen to zero point oscillations 
by the tight transverse trapping. These experimental achievements have revived interest in the theoretical study of 
the properties of ID Bose gases. In most applications, a single parameter, the effective ID scattering length anj, 
is sufficient to describe the interatomic potential, which in this case can be conveniently modeled by a iS-function 
pseudopotential. For repulsive effective interactions the relevant model is provided by the Lieb-Liniger Hamiltonian 
Q . Many properties of this integrable model such as the ground-state energy , the excitation spectrum Q and the 
thermodynamic functions at finite temperature |^ were obtained exactly in the 60' using the Bethc ansatz method. 
Less is known about correlation functions, for which analytic results were obtained only in the strongly interacting 
regime of impenetrable bosons |lCll| and for the long-range behavior of the one-body density matrix . More recently, 
the properties of correlation functions of the Lieb-Liniger model have attracted considerable attention and the short- 
range expansion of the one-body density matrix [T^ , as well as the local two- and three-body correlation function [T^ 
have been investigated. However, a precise determination of the spatial variation of correlation functions for arbitrary 
interaction strength is lacking. 

We use exact quantum Monte Carlo methods to investigate the behavior of correlation functions in the ground 
state of the Lieb-Liniger model. Over a wide range of values for the interaction strength, we calculate the one- and 
two-body correlation function and their Fourier transform giving, respectively, the momentum distribution and the 
static structure factor of the system. These results were already presented in a previous study (Ref. 0|) and are 
briefly reviewed here. We present new results on the local three-body correlation function ranging from the weakly- 
up to the strongly-interacting regime. We have also investigated in more details the long- and short-range behavior 
of the one-body density matrix, including a discussion of finite-size effects and of the validity of the trial function 
used for importance sampling. We always provide quantitative comparisons with known analytical results obtained 
in the weak- or strong-correlation regime, or holding at short or large distances. In the case of the local three-body 
correlator we also compare with available experimental results obtained from three-body loss measurements Q . 

The structure of the paper is as follows. In section^]we introduce the definitions of the spatial correlation functions 
and their Fourier transform. In section ITTll wc discuss the Lieb-Liniger model and give a summary of the main known 
results concerning correlation functions in this model. Section IIVI is devoted to a brief description of the quantum 
Monte Carlo method used for the numerical solution of the Schrodinger equation. The optimization of the trial wave 
function used for importance sampling and the effects due to finite-size are discussed here. The results for correlation 
functions are presented in section Finally, in section Ivll we draw our conclusions. 



II. CORRELATION FUNCTIONS 



We use the first quantization definition of correlation functions in terms of the many-body wave function of the 
system ^{zi, zm), where zi, zm denote the coordinates of the A'' particles. We are interested in the regime where 
the most relevant fluctuations are of quantum nature. In the following we will always consider homogeneous systems 
at T = and the ground-state wave function of the system will be denoted by 4'o(^i, zn)- 
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The one-body density matrix 171(2:) describes the spatial correlations in the wave function and is defined as 

N J'^q{zi+ z,...,ZN)^o{zi,...,ZN)dz2...dzN 
n J \'$o{zi,...,ZN)rdzi...dz]s[ 

where n = N/L is the one-dimensional density. The normalization of .91(2;) is chosen in such a way that 51 (0) = 1. 

Another important quantity is the pair distribution function ,92(1 21 — 22 1), which corresponds to the probability of 
finding two particles separated by |zi — Z2h 

N{N -1) J\^o{zi,...,ZN)\^dz3...dzN 

g2{\Zl-Z2\) = 2 TT— — -T— . (2) 

j \^'o{zi,...,ZN)rdzi...dzN 

The normalization is chosen in such a way that at large distances (72(2) goes to 1 — i.e. becomes unity in the 
thermodynamic limit N —^ 00. 

The value at zero distance of the three-body correlation function gives the probability of finding three particles at 
the same position in space 

^ NjN - -2)1 1^0(0, 0, 0, Z4, zm)\^ dzi...dzN , . 

^""^ ' n3 / 1*0(21,.. .,zjv)Pdzi...dz^ ■ ^' 

The knowledge of the density dependence of ,93(0) allows one to estimate the rate of three-body recombinations 
which is of great experimental relevance Q. 

Much useful information can be obtained from the Fourier transform of the above correlation functions. The 
momentum distribution n{k) is related to the one-body density matrix [Eq. (Q]: 

n{k)=nje''''g,{z)dz, (4) 

and the static structure factor S{k) is instead related to the pair distribution function [Eq. ij^jl]: 

5(fc) = 1 + n /e*'=^(52(2)- l)d2 . (5) 



The momentum distribution can be measured in time of flight experiments and the static structure factor using 
Bragg spectroscopy 0. 

III. LIEB-LINIGER HAMILTONIAN 

A bosonic gas at T = 0, confined in a waveguide or in a very elongated harmonic trap, can be described in 
terms of a one-dimensional model if the energy of the motion in the longitudinal direction is insufficient to excite 
the levels of the transverse confinement. Further, if the range of the interatomic potential is much smaller than the 
interparticle distance and the characteristic length of the external confinement, a single parameter is sufficient to 
describe interactions, namely the effective one-dimensional scattering length aiD- In this case the particle-particle 
interactions can be safely modeled by a (5-pseudopotential. Such a system is described by the Lieb-Liniger (LL) 
Hamiltonian 01 : 

= - ^ II 5^2 + .91^ " ' ('^) 

i—l * i<j 

where the coupling constant gio is related to aio by guj — —2h'^/maio, rn being the particle mass. In the presence 
of a tight harmonic transverse confinement, characterized by the oscillator length a±_ = ^Jh/muo^_^ the scattering 
length was shown to exhibit a non trivial behavior in terms of the 3D s-wave scattering length a^o due to virtual 
excitations of transverse oscillator levels 

aiD = -a± f — - 1.0326 ) . (7) 
In typical experimental conditions, far from a Feshbach resonance, one has 03^1 <Si a±. In this case the above 



equation simplifies to auj = —a^^/a^o, which coincides with the mean-field prediction '|l7j. In the vicinity of a 
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FIG. 1: Energy per particle: Bethe ansatz solution (solid line); GP limit (dashed line); TG limit (dotted line). The circles are 
the results of DMC calculations. Energies are in units of h"^ / {ma\^). 



magnetic Fcshbach resonance the value of a^o can become comparable with a±^ and the coupling constant gio varies 
over a wide range as it goes through a confinement induced resonance . 

All properties in the model depend only on one parameter, the dimensionless density n|aiD|. Contrary to the 3D 
case, where at low density the gas is weakly interacting, in ID small values of the gas parameter ?7|ai£)| correspond to 
strongly correlated systems. This peculiarity of ID systems can be easily understood by comparing the characteristic 
kinetic energy 7l^n^/^/2m, where D denotes the dimensionality, to the mean-field interaction energy gn. In 3D, 
g^on (X n|a3£)| <g; n^/^ if n|a3£)|^ <Si 1. In ID, instead, gion oc n/|ai£)| ^ if n|ai£)| ^ 1. 

The ground-state energy of the Hamiltonian ^ with gio > was first obtained by Lieb and Liniger 01 using the 
Bethe ansatz method. The energy of the system is conveniently expressed as E/N = e{n\ai]j\)h^rv^ /2m, where the 
function e(ri|ai£)|) is obtained by solving a system of integral equations. In the mean-field Gross-Pitaevskii (GP) 
regime, ri|ai£)| ^ 1, the energy per particle is linear in the density E^^ /N ~ gion/2, while in the strongly correlated 
Tonks-Girardeau (TG) regime. n|ai£)| <C 1, the dependence is quadratic E'^'^ /N = /6m. The equation of state 

resulting from a numerical solution of the LL integral equations is shown in Fig. ^ as a function of the gas parameter 
n\aiD\- 

In the TG regime, the energy of an incident particle is not sufficient to tunnel through the repulsive interaction 
potential and two particles will never be at the same position in space. This constraint, together with the spatial 
peculiarity of ID systems, acts as an effective Fermi exclusion principle. Indeed, in this limit, the system of bosons 
acquires many Fermi-like properties. There exists a direct mapping of the wave function of strongly interacting bosons 
onto a wave function of non interacting spinless fermions due to Girardeau The chemical potential, /j, = dE/dN, 
of a TG gas equals the ID Fermi energy = 7i^fc|,/2m, where kp = nn is the Fermi wave vector, and the speed of 
sound c, defined through the inverse compressibility mc^ = nd^/dn, is given by c = Tikp /m. 

Further, the pair distribution function 
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g2{z) = l-- ^, (8) 

which exhibits Friedel-like oscillations, and the static structure factor of the TG gas 

can be calculated exactly exploiting the Bose- Fermi mapping [T^ . 

The one-body density matrix 51(2;) of a TG gas has been calculated in terms of series expansions holding at small 
and large distances in Ref. 0| . The leading long-range term decays as 

gi{z) = ^ = , (10) 



{A = 1.28... is Glaisher's constant) and yields an infrared divergence in the momentum distribution n{k) oc 1/ yj\k\/n. 
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Outside the TG regime full expressions of the correlation functions are not known. The long-range asymptotics 
can be calculated using the hydrodynamic theory of the low-energy phonon-like excitations [Tllll3. l20l | . For gi{z) one 
finds the power-law decay 

^ (11) 

where a = mc/{2TrhTi) and Casympt is a numerical coefficient. This result holds for distances \z\ 3> where ^ = 
h/ (v^TOc) is the healing length. In the TG regime c = Trhn/ni and a = 1/2 as anticipated above. In the opposite GP 
regime (n|aiD| ^ 1) one finds a = l/(TTy^2n\aiD\), yielding a vanishingly small value for a. The power-law decay 
of the one-body density matrix excludes the existence of Bose-Einstein condensation in infinite systems The 
behavior of the momentum distribution for |fc| <^ 1/f follows immediately from Eq. Hll|l : 



Tl{k^ — Casympt 



2n ^- 0Fr(l-f) 

r(f) 



(12) 



where T{z) is the Gamma function. 

The hydrodynamic theory allows one to calculate also the static structure factor in the long-wavelength regime, 
|fc| <§C 1/^. One finds the well-known Feynmann result j23 

Recently, the short-range behavior of the one-, two-, and three-body correlation functions has also been investigated. 
The value at z = of the pair correlation function for arbitrary densities can be obtained from the equation of state 
through the Hcllmann-Feynman theorem [l3| : 

..(0) = -^^^e', (14) 

where the derivative of the equation of state e(n|ai£i|) is to be taken with respect to n|ai£i|. 

The z = value of the three-body correlation function was obtained within a perturbation scheme in the regions 
of strong and weak interactions It is very small in the TG limit (n|aii)| <C 1) 

.93(0) = — , (15) 

and goes to unity in the GP regime (njaiul 3> 1) 

53(0) = 1 - ■ (16) 

The first terms of the short-range expansion of gi{z) can also be calculated from the knowledge of the equation of 
state [III 

.gi(z) 1 - i(e + e'n\a,D\)inz)^ + |("NI') , (17) 
holding for arbitrary densities and for small distances n\z\ <^ 1. 

IV. QUANTUM MONTE CARLO METHOD 

We use Variational Monte Carlo (VMC) and Diffusion Monte Carlo (DMC) methods in order to study the ground- 
state properties of the system. In VMC one calculates the expectation value of the Hamiltonian over a trial wave 
function ■(/'^(R, A, i3, ...), where R = (zi,...,Z7v) denotes the particle coordinates and A, B, ... are variational 
parameters. According to the variational principle, the energy 
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provides an upper bound to the ground-state energy. Evmc ^ Eq. The variational parameters B, ... are optimized 
by minimization of the variational energy (|18|l . 

In order to remove the bias in the estimate of the ground-state energy caused by the particular choice of the trial 
wave function, we resort to the DMC method, which allows one to solve exactly, apart from statistical uncertainty, 
the many-body Schrodinger equation of a Bose system at zero temperature |23l |. The evolution in imaginary time, 
T = it/Ti, is performed for the product /(R, r) = ?/'t(R)5'(R, t), where \I'(R, r) denotes the wave function of the 
system and tPtCR-) is a trial function used for importance sampling. The time-dependent Schrodinger equation for 
the function /(R, r) can be written as 

- ^-^g^'^^ = - D\/lf{R,T) + DVr[F(R)/(R,t)] + [El{II) - i?.e/].f (R, r) , (19) 

where i?i(R) = -0t(R)~^^^V't(R-) denotes the local energy, F(R) = 2-0t(R)~^Vr^/'t(R) is the quantum drift force, 
D ~ fi^ /{2m) plays the role of an effective diffusion constant, and E^ef is a reference energy introduced to stabilize 
the numerical evaluation. The energy and other observables of the state of the system are calculated from averages 
over the asymtpotic distribution function /(R, t — > oo). It is easy to check by decomposing /(R, t) on the basis 
of stationary states of the Hamiltonian, that contributions of the excited states vanish exponentially fast with the 
imaginary time r and asymptotically one obtains limi-^oo /(R, ''') = '0t(R)^o(R') for sl\ trial wave functions non 
orthogonal to the ground-state wave function 'I'o(R). The local energy El sampled over the asymptotic distribution 
equals exactly the ground-state energy 

^.^/i..(R)/(R,r)dR^(^.|^l^o),^^. (20) 
r-oo //(R,r)dR (Vt|*o) ^ ' 

In the present study the trial wave function is chosen of the the Bijl-Jastrow form: 

^T{zu...,ZN)=X{f{z^J), (21) 

where f{z) is a two-body term chosen as 

j Acos[k{\z\-B)], \z\<Z 
^^^>~\s\nP{^\z\/L), \z\>Z. ^^^> 

We consider A'' particles in a box of size L with periodic boundary conditions. In the construction of the trial wave 
function we have ensured that f{z) is uncorrelated at the box boundaries, f{z = ztL/2) = 1 and the derivative 
f'{z = ztL/2) = 0. For \z\ < Z, the Bijl-Jastrow term /(z) corresponds to the exact solution of the two-body problem 
with the potential giD^iz) and provides a correct description of short-range correlations. Long-range correlations 
arising from phonon excitations are instead accounted for by the functional dependence of f{z) for \z\ > Z [T^ . 
The z = boundary condition /'(0+) — /'(0~) = 2/(0)/|ai£)|, which accounts for the (5-function potential, fixes the 
parameter k through the relation fc|aiD| tanfci? = 1. The remaining parameters A, B and /3 are fixed by the continuity 
conditions at the matching point z = Z oi the function f{z), its derivative f'{z) and the local energy — 2/"(z)//(z). 
The value of the matching point Z is a variational parameter which we optimize using VMC. The TG wave function of 
Ref. = rii<j I sin[7r(z,; — Zj)/L]\, is obtained as a special case of our trial wave function ^pT tor Z = B = L/2 

and kL = tt. 

The choice of a good trial wave-function is crucial for the efficiency of the calculation. In order to prove that our 
trial wave-function is indeed very close to the true ground-state ^'o (R)jWe compare in Table I the variational energy 
Evmc with the exact solution based on the use of the Bethe ansatz Q- The corresponding results obtained using 
DMC coincide, within statistical uncertainty, with the exact ones and are shown in Fig. ^ 

Besides the ground-state energy Eq, the DMC method gives exact results also for local correlation functions, such as 
the pair distribution function g2{z) and the three-body correlator 53(0), for which one can use the method of "pure" 
estimators Instead, in the calculation of the non-local one-body density matrix gi{z), the bias from the trial 

wave function can be reduced using the extrapolation technique: (vlfol Ajvl'o) = 2 (4'o|A|'(/;t) — (V'tI^IV't), written here 
for a generic operator A. The "mixed" estimator ("I/ol Aj'i/'T) is the direct output of the DMC calculation and the 
variational estimator (i/jtIAIiPt) is obtained from a VMC calculation. This procedure is accurate only if ipT — ^'o- In 
the present study DMC and VMC give results for gi{z) which are very close and we believe that the extrapolation 
technique removes completely the bias from ^'^(R). 

Calculations are carried out for a finite number of particles N. In order to extrapolate to infinite systems, we 
increase N and study the convergence in the quantities of interest. The dependence on the number of particles 
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n\aiD\ 


Ell/N 


EvMc/N 


10"^ 


1.6408 ■ 10"" 


1.6415(1) ■ lO"*^ 


0.03 


1.3949 ■ 10"=* 


1.3957(3) ■ 10"^ 


0.3 


9.0595 ■ 10"^ 


9.093(1) • 10"^ 


1 


0.5252 


0.5259(1) 


30 


26.842 


27.19(5) 


10^ 


981.15 


983.6(3) 



TABLE I: Energy per particle for different values of the gas parameter ?i|aii3|: Ell — exact result obtained by solving the 
Lieb-Liniger equations, Evmc — variational result, Eq. 1181 . obtained from optimization of the trial wave function 1221 . 
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FIG. 2: Study of finite-size effects in the calculation of the one-body density matrix gi{z) at n|ai_D| = 30. The thin solid line 
corresponds to the power-law decay yl/(n|z|'' '^'') with A obtained from a best fit to the A'' = 500 result. 



(so called finite-size effect) is more pronounced in the regime njaiDj ^ 1, where correlations extend to very large 
distances. Finite size effects can be best investigated by considering the one-body density matrix. In Fig.|21we show 
gi{z) at a fixed density, n|ai£)| = 30, for systems of = 50,100,200,500 particles and we compare the long-range 
behavior with the power-law decay gi{z) oc l/(n|2;|)". For all values of A^ deviations from a power-law decay are 
visible close to the boundary of the box, z = i/2, due to the use of periodic boundary conditions. We find that 
the slope of gi{z) at large distances depends on N , approaching the predicted hydrodynamic value as N increases. 
For A^ ~ 500 wc recover the result a = 0.04. For smaller values of n|aiD|, finite-size effects are less visible and, for 
practical purposes, calculations with A^ < 500 are sufficient. 



V. RESULTS 



We calculate the pair distribution function 52(2) for densities ranging from very small values of the gas parameter 
n|ai£)| ^ 1 (TG regime) up to n|aiD| ^ 1 (GP regime). The results are presented in Fig. 13 In the GP regime 
the correlations between particles are weak and g2{z) is always close to the asymptotic value 52(|-z| 00) ~ 1. By 
decreasing n|ai£)| (thus making the coupling constant gio larger) we enhance beyond-mean field effects and the role 
of correlations. For the smallest considered value of the gas parameter n|ai£)| = 10"'^, we see oscillations in the pair 
distribution function, which is a signature of strong correlations present in the gas. At the same density we compare 
the pair distribution function with the one corresponding to a TG gas, Eq. finding no visible difference. 

In the same figure we show the analytical predictions for the value of ^2(0), Eq. (|14|l . In the TG regime particles 
are never at the same position and consequently (72(0) = 0. With a weaker interaction between particles, we find a 
finite probability that two particles come close to each other in agreement with Eq. (|14|) . As wc go further towards 
the GP regime, the interaction potential becomes more and more transparent and we approach the ideal gas limit 
52(0) = 1. 

In Fig. 21 we present results of the static structure factor obtained from g2{z) according to Eq. 0. At the smallest 
density, njaiul = 10"^, our results are indistinguishable from the S{k) of a TG gas [Eq. ©]. For all densities the small 
wave vector part of S{k) is dominated by phononic excitations. We compare the DMC results with the asymptotic 
linear slope [Eq. 113() ]. We see that in the strongly correlated regime the phononic contribution provides a correct 
description of S(k) up to values of k of the order of the inverse mean interparticle distance n. In the GP regime the 



FIG. 3: Pair distribution function for different values of the gas parameter. In ascending order of the value at zero n|ai_D| = 
10^^, 0.3, 1, 30, lO'*. Arrows indicate the value of 32(0) as obtained from Eq. 11411 . 
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FIG. 4; Static structure factor at the density 7i|ai_D| = 10 0.3, 1, 30, 10"^ (solid lines). The dashed lines are the corresponding 
long- wavelength asymptotics from Eq. (I13II . 

healing length becomes significantly larger than the mean interparticle distance, leading to deviations from the linear 
slope for smaller values of k. 

We show in Fig.[Slthe results for the value at zero distance of the three-body correlation function, Eq. @, calculated 
over a large range of densities. At large density, ri,|aiD| = 10"', the probability of three-body collisions is large and 
the result of Bogoliubov theory, Eq. lfTB|) . provides a good description of 53(0). By reducing n|ai_D| the value of 33(0) 
decreases, becoming vanishingly small for values of the gas parameter n|ai£)| ^ 1. In order to resolve the dependence 
of (73(0) on the density in the TG regime, we plot the results on a log- log scale (inset of Fig. ISJ. We obtain that 
53(0) is proportional to the fourth power of the gas parameter in agreement with Eq. (|15|l . A reliable evaluation of 
the three-body correlator for small densities is difficult due to the very small value of 53 (0) itself. It is interesting to 
notice that gfCO) is close to 33(0) over the whole density range. This estimate of (73(0) has been discussed in Ref. 0]. 
The coefficient of three-body losses has been measured in quasi-lD configurations realized with deep two-dimensional 
optical lattices Q . The value of (73 (0) extracted from these measurements is also shown in Fig. |S| 

We calculate the spatial dependence of the one-body density matrix, Eq.^ for different values of the gas parameter. 
At small distances we compare the DMC results with the short-range expansion, Eq. 117|) . finding agreement for nz 
(see Fig.E)). For distances larger than the healing length ^ we expect the hydrodynamic theory to provide a correct 
description. The long-range decay shown in Fig. |7| exhibits a power-law behavior in agreement with the prediction of 
Eq. (|11|) . The coefficient of proportionality in Eq. is fixed by a best fit to the DMC results. The small deviations 
from the power law-decay at the largest distances [z k, L/2) are due to finite size effects (see Fig.|2l. 

In the weakly interacting GP regime the coefficient Casympt of Eq. Hll|l can be calculated from a hydrodynamic 



approach One obtains 




(23) 



where 7 = 0.577 is Euler's constant and a ~ 'mc/{2Trhn). 
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FIG. 5: Value at zero distance of the three-body correlation function g3{0) (circles); GP limit, Eq. 11611 . dashed line; mean-field 
factorization, (73(0) = ('72(0))'' solid line. Inset: small density region on a log-log scale: TG limit, Eq. 1151 . dotted line. Open 
symbol: measured value of (73(0) from three-body loss experiments U 





pDMC 
^ asympt 


asympt 


C asympt 


1000 


1.02 


1.0226 


1.0226 


30 


1.06 


1.0588 


1.0579 


1 


0.951 


0.9646 


0.9480 


0.3 


0.760 


0.8145 


0.7814 


0.001 


0.530 


0.5746 


0.5227 



TABLE II: Coefficient of the long-range decay of the one-body density matrix as defined in Eq. IIH . The first column is the 
one- dimensional gas parameter, the second column is the coefficient extracted from the best fit to the DMC results (see Fig. [3, 
the third column is Popov's prediction, and the fourth column is Eq. 12311 . The value of the gas parameter njaiijl = 0.001 
corresponds to the deep TG regime where one can apply Eq. 1101 yielding Ca^^mpt ~ 0.5214. 




Although result H23|l is formally derived in the weakly interacting Umit, a ^ 1, it works well in the whole range 
of densities. Indeed, the value of Casympt obtained from the best fit to gi{z), as shown in Fig.[71 is always in good 
agreement with the prediction (|23|l . For example, in the strongly-interacting TG regime the comparison between Eg . 
(|23|l and the exact result in Eq. 1)11) |l gives only 0.3% difference. A different expression was obtained by Popov |2g 

(and later recovered in j23|) giving c£,°^°p( = • Both expressions coincide for small values of a, but Popov's 

coefficient yields larger errors approaching the TG regime, with 10% maximal error, as it was pointed out in Ref. |28| . 
A comparison of the different coefficients is presented in Table II. 

The momentum distribution, Eq. Q), is obtained from the Fourier transform of the calculated one-body density 




FIG. 6: Short range behavior of the one-body density matrix at density n\aiD \ = 10 ^, 0.3, 1, 30 (solid lines), compared with 
the series expansion at small distances [Eq. 1171 1 (dashed lines). 
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FIG. 7: Long-range behavior of the one-body density matrix (sohd hnes), best fits to the long- wavelength asymptotics from 
Eq. mn (dashed lines). Values of the density are n\a\D I = 10"^0.3, 1,30, 10^ The arrows indicate the value of the product 
of the density and the healing length ^n: the leftmost corresponds to n|aiD| = 10~'^, the rightmost to n|ai_D| ~ 10^ 
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FIG. 8: Momentum distribution at density n|ai_D| = 10^"^, 0.3, 1, 30, 10''. The dashed lines correspond to the infrared behavior 
of Eq. 11121 . The arrows indicate the value of l/^n: the rightmost corresponds to n|ai_D| = 10^^, the leftmost to n|aii5| — 10''. 

matrix at short distances and the best fitted power-la"w decay at large distances. The momentum distribution exhibits 
the infrared divergence of Eq. H12|) . We present the results for n{k) by plotting in Fig.|Slthe combination kn{k), where 
the divergence is absent. We notice that the infrared asymptotic behavior is recovered for values of k considerably 
smaller than the inverse healing length 1/^. At large k the numerical noise of our results is too large to extract 
evidences of the 1 /fc^ behavior predicted in Ref . . 

VI. CONCLUSIONS 

This paper presents a thorough study of the correlation functions in a one-dimensional homogeneous Bose gas 
described by the Lieb-Liniger Hamiltonian. The correlation functions are calculated for all interaction regimes using 
exact quantum Monte Carlo methods. The results on the pair distribution function, one-body density matrix and 
their Fourier transformations have already been presented in Ref. |l4j and are briefly revie^wed here. We carry out a 
more detailed study of the short- and long-range behavior of the one-body density matrix, including the comparison 
with analytic expansions. We also calculate the probability of finding three particles at the same spatial position 
as a function of the interaction strength and we compare it with the asymptotic results holding in the weakly- and 
strongly-interacting regime and with available experimental results from three-body loss measurements. 
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